//
//  modules.c
//  PyroTools
//
//  Created by Timothée Poisot on 12/07/11.
//  Copyright 2011 Université Montpellier 2. All rights reserved.
//
//  General functions for analysis
//

#include <stdio.h>
#include <math.h> 
#include <stdlib.h>
#include <string.h>
#include <time.h>

#include "modules.h"
#include "stat-functions.h"
#include "genet-functions.h"




int newoutfile (const char * argv[])
{
    FILE *outfile;
    outfile = fopen(argv[2], "w");
    printf("New output file %s created\n\n",argv[2]);
    fclose(outfile);
    return 0;
}

int betadiv (const char * argv[])
{
    FILE *pop1, *pop2, *seq1, *seq2, *output;
    char str[MaxPop];
    

    int size1 = 0;
    int size2 = 0;
        
    pop1=fopen(argv[2], "r");
    pop2=fopen(argv[3], "r");
    
    // START number of OTUs
    
    if (pop1) {
        while (fscanf(pop1, "%i", str)!=EOF)
            ++size1;
    }
    if (pop2) {
        while (fscanf(pop2, "%i", str)!=EOF)
            ++size2;
    }
    
    fclose(pop1);
    fclose(pop2);
    
    printf("Sample 1:\t%i unique OTUs\nSample 2:\t%i unique OTUs\n\n",size1,size2);   
    
    //
    
    // ARRAYS to hold population size and distances
    int* ps1 = (int*)malloc(size1 * sizeof(int*));
    int* ps2 = (int*)malloc(size2 * sizeof(int*));
    double* dist = (double*)malloc(size1 * size2 * sizeof(double*));
    double* weight = (double*)malloc(size1 * size2 * sizeof(double*));
    
    pop1=fopen(argv[2], "r");
    pop2=fopen(argv[3], "r");
    
    int popPos1 = 0;
    int popPos2 = 0;
    
    char line[MaxPop];
    int limit;
    
    while (fgets(line,MaxPop,pop1)!=NULL)
    {
        
        sscanf(line,"%i", &limit);
        ps1[popPos1] = limit;
        popPos1++;
        
    }
    
    while (fgets(line,MaxPop,pop2)!=NULL)
    {
        
        sscanf(line,"%i", &limit);
        ps2[popPos2] = limit;
        popPos2++;
        
    }
    
    fclose(pop1);
    fclose(pop2);
    
    int TP1 = iSum(ps1,size1);
    int TP2 = iSum(ps2,size2);
    
    double fs1[size1];
    double fs2[size2];
    
    for(int i = 0;i<size1;++i)
    {
        fs1[i] = (double) ps1[i]/TP1;
    }
    for(int i = 0;i<size2;++i)
    {
        fs2[i] = (double) ps2[i]/TP2;
    }
    
    
    
    
    
    /////////////////////////////
    // START genetic distances
    
    // Read sequences 1
    seq1=fopen(argv[4], "r");
    // Dynamic allocation of array of sequences 1
    char** se1;
    char* temp1;
    
    se1 = (char**)malloc(size1 * sizeof(char*));
    temp1 = (char*)malloc(size1 * MaxSL * sizeof(char));
    for (int i = 0; i < size1; i++) {
        se1[i] = temp1 + (i * MaxSL);
    }
    //
    
    int i1 = 0;
    if ( seq1 != NULL )
    {
        char l1 [ MaxSL ]; /* or other suitable maximum line size */
        
        while ( fgets ( l1, MaxSL, seq1 ) != NULL ) /* read a line */
        {
            for(int i=0;i<MaxSL;++i)
            {
                se1[i1][i] = l1[i];
            }
            ++i1;
        }
        fclose ( seq1 );
    }
    
    // Read sequences 2
    
    seq2=fopen(argv[5], "r");
    // Dynamic allocation of array of sequences 2
    char** se2;
    char* temp2;
    
    se2 = (char**)malloc(size2 * sizeof(char*));
    temp2 = (char*)malloc(size2 * MaxSL * sizeof(char));
    for (int i = 0; i < size2; i++) {
        se2[i] = temp2 + (i * MaxSL);
    }
    //
    
    int i2 = 0;
    if ( seq2 != NULL )
    {
        char l2 [ MaxSL ];
        
        while ( fgets ( l2, MaxSL, seq2 ) != NULL )
        {
            for(int i=0;i<MaxSL;++i)
            {
                se2[i2][i] = l2[i];
            }
            ++i2;
        }
        fclose ( seq2 );
    }
    
//    output = fopen(argv[6],"w");
    
    int currentSeqPos = 0;
    
    time_t now;
    time(&now);
    printf("Execution started on %s", ctime(&now));
    
    clock_t start = clock();
    for(int pos1 = 0;pos1<size1;++pos1)
    {
        for(int pos2 = 0;pos2<size2;++pos2)
        {
            double di = genetdist(se1[pos1], se2[pos2], MaxSL);
            dist[currentSeqPos] = di;
            weight[currentSeqPos] = fs1[pos1] + fs2[pos2];
            ++currentSeqPos;
        }
        if(pos1==3)
        {
            printf("Estimated runtime:\t%f seconds\n", ((((double)clock() - start) / CLOCKS_PER_SEC) * (size1-3))/4);
        }
    }
    
    // Clock
    time_t endnow;
    time(&endnow);
    printf("Completion runtime:\t%f seconds\n", (((double)clock() - start) / CLOCKS_PER_SEC));
    printf("Execution ended on %s", ctime(&endnow));
    
    // Generate and compute results
    printf("\n-------\nResults\n\n");
    double dme = dMean(dist, currentSeqPos);
    printf("Mean             %g\n",dme);
    double dwme = dWMean(dist, weight,currentSeqPos);
    printf("Weighted mean    %g\n\n",dwme);
    double vari = dVar(dist, currentSeqPos);
    printf("Variance         %g\n",vari);
    double wvari = dWSd(dist, weight, currentSeqPos);
    printf("W. Variance      %g\n\n",wvari);
    double cov = dCoeffVar(dist, currentSeqPos);
    printf("Var. coeff       %g\n",cov);
    double wcov = dWCV(dist, weight, currentSeqPos);
    printf("W. Var. coeff    %g\n",wcov);
    printf("-------\nResults written to file\n");
    output = fopen(argv[6],"a");
    fprintf(output,"%s %s %g %g %g %g %g %g\n",argv[2],argv[3],dme,dwme,vari,wvari,cov,wcov);
    fclose(output);
    printf("-------\n\n\n");
    
    return 0;
}